      Subroutine RK4_Ghost
      include "param.fi"
      include "common.fi"

      Integer I, ILevelRK4

      Double Precision kRho(NTOT3T, 4), oRho(NTOT3T)
      Double Precision kU1(NTOT3T, 4), oU1(NTOT3T)
      Double Precision kU2(NTOT3T, 4), oU2(NTOT3T)
      Double Precision kU3(NTOT3T, 4), oU3(NTOT3T)
      Double Precision kX1(NTOT3T, 4), oX1(NTOT3T)
      Double Precision kX2(NTOT3T, 4), oX2(NTOT3T)
      Double Precision kX3(NTOT3T, 4), oX3(NTOT3T)
      Double Precision dTempStep(4), dCoefStep(4)

      Do I = 1, NTOT3T
        oRho(I) = Rho(I)
        oU1(I) = U1(I)
        oU2(I) = U2(I)
        oU3(I) = U3(I)
        oX1(I) = X1(I)
        oX2(I) = X2(I)
        oX3(I) = X3(I)
      Enddo

      dTempStep(1) = 0.0
      dTempStep(2) = 0.5
      dTempStep(3) = 0.5
      dTempStep(4) = 1.0
      dCoefStep(1) = 1.0 / 6.0
      dCoefStep(2) = 1.0 / 3.0
      dCoefStep(3) = 1.0 / 3.0
      dCoefStep(4) = 1.0 / 6.0

      Do ILevelRK4 = 1, 4

      If (ILevelRK4 .GT. 1) Then
      Do I = 1, NTOT3T
        Rho(I) = oRho(I) 
     &         + dTempStep(ILevelRK4) * DT * kRho(I, ILevelRK4 - 1)
        U1(I) = oU1(I) 
     &         + dTempStep(ILevelRK4) * DT * kU1(I, ILevelRK4 - 1)
        U2(I) = oU2(I) 
     &         + dTempStep(ILevelRK4) * DT * kU2(I, ILevelRK4 - 1)
        U3(I) = oU3(I) 
     &         + dTempStep(ILevelRK4) * DT * kU3(I, ILevelRK4 - 1)
      If (ALE .EQ. ALELagrange) Then
        X1(I) = oX1(I)
     &         + dTempStep(ILevelRK4) * DT * kX1(I, ILevelRK4 - 1)
        If (X1(I).GT.XMAX) X1(I) = X1(I) - CUBESIZE
        If (X1(I).LT.XMIN) X1(I) = X1(I) + CUBESIZE
        If (DIMM.EQ.DIMM3D) Then
          X2(I) = oX2(I)
     &         + dTempStep(ILevelRK4) * DT * kX2(I, ILevelRK4 - 1)
          If (X2(I).GT.XMAX) X2(I) = X2(I) - CUBESIZE
          If (X2(I).LT.XMIN) X2(I) = X2(I) + CUBESIZE
        Endif
        X3(I) = oX3(I) 
     &         + dTempStep(ILevelRK4) * DT * kX3(I, ILevelRK4 - 1)
        If (X3(I).GT.XMAX) X3(I) = X3(I) - CUBESIZE
        If (X3(I).LT.XMIN) X3(I) = X3(I) + CUBESIZE
        If (tagRegroup .NE. 0) Then
          Call Regroup
        Endif
      Endif
      Enddo
      Endif


      If (Shepard .EQ. 1) Then
        Call Calc_ShepardCoef
      Endif

      Call PressureField

      Call Calc_dRho_dT
      Do I = 1, NTOT3T
        kRho(I, ILevelRK4) = Temp1(I)
      Enddo

      Call Calc_dU_dT
      Do I = 1, NTOT3T
        kU1(I, ILevelRK4) = Temp1(I)
        If (DIMM.EQ.DIMM3D) kU2(I, ILevelRK4) = Temp2(I)
        kU3(I, ILevelRK4) = Temp3(I) + Gravity
      Enddo

      If (ALE .EQ. ALEEuler) Then
        Call Grad(U1, DU1DX, DU1DY, DU1DZ)
        If (DIMM.EQ.DIMM3D) Call Grad(U2, DU2DX, DU2DY, DU2DZ)
        Call Grad(U3, DU3DX, DU3DY, DU3DZ)
        Call Calc_NLTerm
        Do I = 1, NTOT3T
          kU1(I, ILevelRK4) = kU1(I, ILevelRK4) + Temp1(I)
          If (DIMM.EQ.DIMM3D) kU2(I, ILevelRK4) =
     &                        kU2(I, ILevelRK4) + Temp2(I)
          kU3(I, ILevelRK4) = kU3(I, ILevelRK4) + Temp3(I)
        Enddo
      Endif

      Call Calc_VisTerm
      Do I = 1, NTOT3T
        kU1(I, ILevelRK4) = kU1(I, ILevelRK4) + Temp1(I)
        If (DIMM.EQ.DIMM3D) kU2(I, ILevelRK4) =
     &                      kU2(I, ILevelRK4) + Temp2(I)
        kU3(I, ILevelRK4) = kU3(I, ILevelRK4) + Temp3(I)
      Enddo

      If (ALE .EQ. ALELagrange) Then
        Do I = 1, NTOT3T
          kX1(I, ILevelRK4) = U1(I) 
     &       + dTempStep(ILevelRK4) * DT * kU1(I, ILevelRK4)
          kX2(I, ILevelRK4) = U2(I)
     &       + dTempStep(ILevelRK4) * DT * kU2(I, ILevelRK4)
          kX3(I, ILevelRK4) = U3(I)
     &       + dTempStep(ILevelRK4) * DT * kU3(I, ILevelRK4)
        Enddo
      Endif

      Enddo ! ILevelRK4

      Do I = 1, NTOT3T
        Rho(I) = oRho(I)
     &         + DT / 6.0 * kRho(I, 1)
     &         + DT / 3.0 * kRho(I, 2)
     &         + DT / 3.0 * kRho(I, 3)
     &         + DT / 6.0 * kRho(I, 4)
        U1(I) = oU1(I) 
     &         + DT / 6.0 * kU1(I, 1)
     &         + DT / 3.0 * kU1(I, 2)
     &         + DT / 3.0 * kU1(I, 3)
     &         + DT / 6.0 * kU1(I, 4)
        U2(I) = oU2(I) 
     &         + DT / 6.0 * kU2(I, 1)
     &         + DT / 3.0 * kU2(I, 2)
     &         + DT / 3.0 * kU2(I, 3)
     &         + DT / 6.0 * kU2(I, 4)
        U3(I) = oU3(I) 
     &         + DT / 6.0 * kU3(I, 1)
     &         + DT / 3.0 * kU3(I, 2)
     &         + DT / 3.0 * kU3(I, 3)
     &         + DT / 6.0 * kU3(I, 4)
      If (ALE .EQ. ALELagrange) Then
        X1(I) = oX1(I) 
     &         + DT / 6.0 * kX1(I, 1)
     &         + DT / 3.0 * kX1(I, 2)
     &         + DT / 3.0 * kX1(I, 3)
     &         + DT / 6.0 * kX1(I, 4)
        If (X1(I).GT.XMAX) X1(I) = X1(I) - CUBESIZE
        If (X1(I).LT.XMIN) X1(I) = X1(I) + CUBESIZE
        If (DIMM.EQ.DIMM3D) Then
          X2(I) = oX2(I) 
     &         + DT / 6.0 * kX2(I, 1)
     &         + DT / 3.0 * kX2(I, 2)
     &         + DT / 3.0 * kX2(I, 3)
     &         + DT / 6.0 * kX2(I, 4)
          If (X2(I).GT.XMAX) X2(I) = X2(I) - CUBESIZE
          If (X2(I).LT.XMIN) X2(I) = X2(I) + CUBESIZE
        Endif
        X3(I) = oX3(I) 
     &         + DT / 6.0 * kX3(I, 1)
     &         + DT / 3.0 * kX3(I, 2)
     &         + DT / 3.0 * kX3(I, 3)
     &         + DT / 6.0 * kX3(I, 4)
        If (X3(I).GT.XMAX) X3(I) = X3(I) - CUBESIZE
        If (X3(I).LT.XMIN) X3(I) = X3(I) + CUBESIZE

        If (tagRegroup .NE. 0) Then
          Call Regroup
        Endif
      Endif
      Enddo

      Call PressureField

      End Subroutine
